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We compute the next-to-leading order QCD predictions for the vertical flux of atmospheric muons 
and neutrinos from decays of charmed particles, for different PDF's (MRS-R1, MRS-R2, CTEQ- 
4M and MRST) and different extrapolations of these at small partonic momentum fraction x. We 
find that the predicted fluxes vary up to almost two orders of magnitude at the largest energies 
studied, depending on the chosen extrapolation of the PDF's. We show that the spectral index 
of the atmospheric leptonic fluxes depends linearly on the slope of the gluon distribution function 
at very small x. This suggests the possibility of obtaining some bounds on this slope in "neutrino 
telescopes", at values of x not reachable at colliders, provided the spectral index of atmospheric 
. leptonic fluxes could be determined. 
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The flux of atmospheric neutrinos and muons at very high energies, above 1 TeV, originates primarily from semilcp- 
tonic decays of charmed particles instead of pions and kaons, which are the dominant decay modes at lower energies 
(see for example jjj). This flux is one of the most important backgrounds for "neutrino telescopes", limiting their 
sensitivity to astrophysical signals, especially for future km 3 detectors which might be able to observe neutrinos and 
muons at extremely high energies, even up to 10 12 GeV. 

We use perturbative QCD (pQCD), the theoretically preferred model, to compute the charm production. We 
perform a true next-to-leading order (NLO) pQCD analysis of the production of charmed particles in the atmosphere, 
together with a full simulation of the particle cascades down to the final muons and neutrinos. This is done by 
combining the NLO pQCD calculations of charm production and computer routines of Mangano, Nason, and Ridolfi 
HI (called MNR in the following) with the computer simulations of the cascades generated by P YTHIA |4| . These are 
the same programs currently used to compare pQCD predictions with experimental data in accelerator experiments. 

We have already presented results of our calculations in a previous paper § (called GGV1 from now on), in which 
all the details of the program we use can be found. The main goal of our first paper was to compare the fluxes obtained 
with the NLO and the leading order (LO) calculations, i.e. we computed the K factor for the neutrino and muon 
fluxes. This was done to improve on the first study of atmospheric fluxes based on pQCD, performed by Thunman, 
Ingelman, and Gondolo a few years ago in Ref. || (called TIG in the following). TIG used the LO charm production 
cross section computed by PYTHIA, multiplied by a constant K factor of 2 to bring it in line with the NLO values, 
and supplemented by parton shower evolution and hadronization according to the Lund model. 

In GGV1 we found the K factors for different parton distribution functions (PDF's), as function of energy, to be 
in a range between 2.1 and 2.5. A similar analysis was recently made in Pasquali, Reno, and Sarcevic (called PRS 
from now on), with results compatible with ours, using a treatment of the problem complementary to ours. In fact, 
PRS used approximate analytic solutions to the cascade equations in the atmosphere, also introduced by TIG, while 
we make instead a full simulation of the cascades. 

In GGV1 we showed that the approach used by TIG (i.e. multiplying the LO fluxes by an overall K factor of 2) 
was essentially correct, except for their relative low K factor (since K values of 2.2 - 2.4, depending slightly on the 
PDF, provide estimates of the NLO within about 10%). However, while TIG found neutrino and muon fluxes lower 
than the lowest previous estimate, we found instead larger fluxes (by factors of 3 to 10 at the highest energies, about 
10 9 GeV), in the bulk part of previous predictions. The main reason for this difference is studied in this paper. 
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Here we explore the dependence of the atmospheric fluxes on the extrapolation of the gluon PDF at very small 
partonic momentum fraction x, x < 10~ 5 , which is crucial for the fluxes at high energies. As explained below, the 
relevant momentum fraction x of the interacting atmospheric parton is of the order of the inverse of the leptonic 
energy Ei (in the atmospheric rest frame) in GeV. This energy, in turn, is of the order of 0.1 E, where E is the energy 
per nucleon of the incoming cosmic ray in the lab. frame (the atmospheric rest frame). Thus, for E\ > 10 5 GeV, we 
need the PDF's at x < 1CP 5 , values of x which are not reached experimentally. The final fluxes depend mostly on 
the gluon PDF, since this is by far the dominant one at these small x values and charm is mostly produced through 
gluon-gluon fusion processes. 

A concern that has been expressed to us several times is the applicability of the MNR NLO-pQCD calculations, 
mostly done for accelerator physics, to the different kinematic domain of cosmic rays. In response we remark that, 
for the less steep extrapolations of the gluon structure function g(x) that we use at small x, we have large logarithms, 
known as "ln(l/x)" terms, where x ~ y^mf/s, s is the hadronic center of mass squared energy and this x is the 
average value of the hadron energy fraction needed to produce the cc pair. With the extrapolation g(x) ~ x x ~ 1 (see 
below) and A close to 0.5, and possibly for the intermediate choices of A also, there should be no large logarithm. 
The problem arises for A too close to zero. We will attempt to deal with this problem in future work. Moreover, 
contrary to the case in accelerators, we do not have the uncertainty present in the differential cross sections j^] when 
kx is much larger than m c , due to the presence of large logarithms of (k^ + m^)/m^. Because we do not have here 
a forward cut in acceptance, the characteristic transverse charm momentum in our simulations is of the order of the 
charm mass, kx — 0(m c ). 

In this paper, as in GGV1, the MNR program is used to compute the inclusive charm cross section and the cascades 
simulated by PYTHIA are initiated by a single c quark. This is the 'single' mode described in our previous paper 
GGV1, where we argued its advantages. We explained there our normalization of the NLO charm production cross 
sections in the MNR program, and described in detail the computer simulations used to calculate the neutrino and 
muon fluxes, which we briefly review in Sections || and [II. Except for the inclusion of the NLO calculations our 
model closely follows TIG. In Section [V we show the neutrino and muon fluxes we obtain for different low x behaviors 
of the gluon PDF and we compare them with the TIG fluxes. In Section [v[ we give analytic arguments that explain 
and support our results. 

Finally, as in GGV1 (and TIG), we consider only vertical showers for simplicity. We intend to study those from all 
directions in the future. 



II. CHARM PRODUCTION IN PQCD AND CHOICE OF PDF'S 



Our NLO calculation is based on the MNR computer code. The NLO cross section for charm production depends on 
the choice of the parton distribution functions and on three parameters: the charm quark mass m c , the renormalization 
scale fiR, and the factorization scale [ip. In order to calibrate the charm production routines we fit the most recent 
experimental data f|-[ll| (differential and total cross sections) with one and the same combination of m c , (ir, and 
fip, for each PDF we use (see Q for complete details). Several choices of m c , [ir and fip may work equally well. In 
fact the cross sections increase by decreasing fj,p, fj,R or m Cl so changes in the three variables can be played against 
each other to obtain practically the same results. We use just one such choice for each PDF. We intend to further 
study the uncertainty related to this range of possible choices in the future. 

As in GGV1, here we use the PDF's MRS Rl, R2 |l| and CTEQ 4M |]|, with the following parameters. We 
choose /in — rriT, \if — 2m,T for all sets, where mj- is the transverse mass, mx = \Jk T + m^, and 

m c = 1.185 GeV for MRS Rl, (1) 
m c = 1.31 GeV for MRS R2, (2) 
m c = 1.27 GeV for CTEQ 4M. (3) 

The data we use for this 'calibration' of the MNR program are shown in Table 1 and Table 2 of GGV1. In this 
paper, we add to our list of PDF's the latest of the MRS set, the MRST [jl4|, with charm mass 

m c = 1.25 GeV for MRST, (4) 

obtained with the same procedure used for the other PDF's. 

As we will see clearly in Sect. [v|, due to the steep decrease with increasing energy of the incoming flux of cosmic 
rays, only the most energetic charm quarks produced count, and these come from the interactions of projectile partons 
carrying a large fraction of the incoming nucleon momentum. Thus, the characteristic x of the projectile parton, that 
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we call Xi, is large. It is x\ ~ 0(10 1 ). We can, then, immediately understand that very small partonic momentum 
fractions are needed in our calculation, because typical partonic center of mass energies \/I are close to the cc threshold, 
2m c ~ 2 GeV (since the differential cross section decreases with increasing s) while the total center of mass energy 
squared is s = 2m^E (with m^r the nucleon mass, ~ 1 GeV). Calling X2 the momentum fraction of the target 
parton (in a nucleus of the atmosphere), then x\x 2 = s/s — 4m 2 /(2m,/v-E) ~ GeV /E. Thus, x 2 — O(GeV/0.1-E), 
where E is the energy per nucleon of the incoming cosmic ray in the lab. frame. The characteristic energy E c of the 
charm quark and the dominant leptonic energy Ei in the fluxes are Ei ~ E c ~ 0.1-E, thus x 2 — 0(GeV / E{). Namely 
x 2 ^ 10~ 6 , 1(T 7 at Ei ~ 1, 10 PeV. 

For x > 10 -5 (E < 10 3 TeV), PDF's are available from global analyses of existing data. We use four sets of PDF's. 
Three of these, MRS Rl, MRS R2 [§Jj and CTEQ 4M || (used also in GGV1), incorporate most of the latest HERA 
data and cover the range of parton momentum fractions x > 10~ 5 and momentum transfers Q 2 > 1.25 — 2.56 GeV 2 . 
MRS Rl and MRS R2 differ only in the value of the strong coupling constant a s at the Z boson mass: in MRS Rl 
a a (M|) = 0.113 , and in MRS R2 a a (M§) = 0.120. The former value is suggested by "deep inelastic scattering" 
experiments, and the latter by LEP measurements. This difference leads to different values of the PDF parameters 
at the reference momentum Q 2 , = 1.25 GeV 2 , where the QCD evolution of the MRS Rl and R2 PDF's is started. 
The CTEQ 4M is the standard choice in the MS scheme in the most recent group of PDF's from the CTEQ group 
(a a (Af§) = 0.116 for CTEQ 4M). In this paper we also use the very recent MRST @. This new PDF set includes 
all the latest experimental measurements that have become available and, for the first time, an investigation of the 
uncertainty in the gluon distribution function. We will use the main choice of the MRST set, the "central gluon" 
MRST, the central value of the gluon PDF's of the package, which is considered the optimum global choice of this new 
set. The range in Q 2 and x of MRST set is the same as for the older MRS R1-R2 (x > 10~ 5 and Q 2 > 1.25 GeV 2 ), 
and a a (M|) = 0.1175. 

For x <C 1, all these PDF's go as 

xfi(x, Q 2 ) ~ AiX- Xi( > Q2 \ (5) 

where i denotes valence quarks u v ,d Vl sea quarks S, or gluons g. The PDF's we used have Xs(Q 2 ) ^ ^g(Qo); m 
contrast to older sets of PDF's which assumed an equality. As x decreases the density of gluons grows rapidly. At 
x ~ 0.3 it is comparable to the quark densities but, as x decreases it increasingly dominates over them. Quark 
densities become negligible at x < 10~ 3 . 

The PDF's need to be extrapolated to x < 10 -5 (E > 10 3 TeV). Extrapolations based on Regge analysis usually 



propose xg(x) ~ x with A ~ 0.08 |15|, while evolution equations used to resum the large logarithms a s ln(l/x) 
mentioned before, such as the BFKL (Balitsky, Fadin, Kuraev, Lipatov |l6|) find also xg(x) ~ x~ x , but with A ~ 0.5. 

In this work we use extrapolations with different values of A. For the older MRS R1-R2 and CTEQ 4M we consider 
only the two extreme behaviors and the intermediate one that we used in GGV1, namely: (i) a constant extrapolation 
X g (Q 2 ) = for x < 10~ 5 ; (ii) a linear extrapolation of In g(x) as a function oflnx, \ng(x) = — (A g (Q 2 ) + l) lnx + ln A g , 
where X g (Q 2 ) is taken at x = 10~ 5 , the smallest x for which the PDF's are provided (we call A(R1), A(R2) or A(4M) 
the A's so obtained); (iii) an extrapolation with X g (Q 2 ) = 0.5 for x < 10~ 5 . Cases (i) and (iii) are extreme choices 
theoretically justified before |l5|], while (ii) is somewhat in between, with a resulting A ~ 0.2 — 0.3. 

For the new MRST we have included several values of A, in order to test the dependence on this parameter in a 
more complete way: (i) extrapolations with different A's, i.e. X g (Q 2 ) = 0,0.1,0.2,0.3,0.4,0.5 for x < 10 -5 ; (ii) we 
also included the linear extrapolation of In g(x) as a function of In a;, similar to the second intermediate choice of the 
previous list; we will call A(T) the A obtained in this way. 



III. SIMULATION OF PARTICLE CASCADES IN THE ATMOSPHERE 



In this section we briefly describe the computer simulation used to calculate the neutrinos and muons fluxes; a 
more detailed description can be found in GGV1 ||. The charm production process in the atmosphere and the 
particle cascades are simulated by modifying and combining together two different programs: the MNR routines 
and PYTHIA 6.115 |]. 

The MNR program was modified to become an event generator for charm production at different heights in the 
atmosphere and for different energies of the incoming primary cosmic rays. 

The charm quarks (and antiquarks) generated by this first stage of the program are then fed into a second part 
which handles quark showering, fragmentation and the interactions and decays of the particles down to the final 
leptons. The cascade evolution is therefore followed throughout the atmosphere: the muon and neutrino fluxes at sea 
level are the final output of the process. 
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In order to make our results comparable to those of TIG, we keep the same modeling of the atmosphere and of the 
primary cosmic ray flux as in TIG and the same treatment of particle interactions and decays in the cascade. 

We recall however that our main improvements are the inclusion of a true NLO contribution for charm production, 
the use of updated PDF's and, in this second paper, the different extrapolations used for the gluon PDF at low x. 

In the rest of this section we review briefly the model for the atmosphere and the primary flux used in this study, 
which is the same of GGV1 and was introduced originally by TIG. 

We assume a simple isothermal model for the atmosphere. Its density at vertical height h is 

with the parameters, scale height ho = 6.4 km and column density X$ = 1300 g/cm 2 at h — 0, chosen as in TIG 
to fit the actual density in the range 3 km < h < 40 km, important for cosmic ray interactions. Along the vertical 
direction, the amount of atmosphere traversed by a particle, the depth X, is related to the height h simply by 

/>oo 

X= p{h')dh' =X a e- h ' ha . (7) 
hi 

The atmospheric composition at the important heights is approximately constant: 78.4% nitrogen, 21.1% oxygen and 
0.5% argon with average atomic number (A) = 14.5. 

Following TIG ||, we neglect the detailed cosmic ray composition and consider all primaries to be nucleons with 
energy spectrum 



(25,0) 



nucleons 



cm 2 s sr GeV /A 



= foE-'- 1 = (8) 

f 1.7(£7GeV)- 2 - 7 for E < 5 10 6 GeV 
\ 17i(E/GeV)- 3 for E > 5 10 6 GeV 



The primary flux is attenuated as it penetrates into the atmosphere by collisions against the air nuclei. An 
approximate expression for the intensity of the primary flux at a depth X is (see || again) 

4> N {E,X) =e~ x ' AN (f> N (E,0) . (9) 

The nuclear attenuation length A^r, defined as 

has a mild energy dependence through Ajv and Znn, the spectrum- weighted moment for nucleon regeneration in 
nucleon-nucleon collisions. We use the Z^n values in Fig. 4 of Ref. ||. The interaction thickness An is 

\ N (E,h) = = P{ ?1 (11) 

L A ^NAiEjriA^h) 

where nji(h) is the number density of air nuclei of atomic weight A at height h and <jna(E) is the total inelastic 
cross section for collisions of a nucleon N with a nucleus A. This cross section scales essentially as A 2 / 3 , a^A^E) = 
A 2 / 3 unn{E). For unn{E) we use the fit to the available data in Ref. |l7|| . Using our height independent atmospheric 
composition, we simplify Eq. ( pd| ) as follows, 

X N (E, h) = -^L- = 2.44 . (12) 

(A 2 ^) (jnn{E) (jnn(E) 

Here ( ) denotes average and u is the atomic mass unit, that we write as 

u = 1660.54 mb g/cm 2 . (13) 

Therefore in our approximations Xn(E) is independent of height. 
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IV. NEUTRINO AND MUON FLUXES 



We present here the results of our simulations with all the PDF's and the values of A described in Section |1 
The NLO total inclusive charm- ant icharm production cross sections a cS for our four different PDF's are shown in 
Fig. 1 over the energy range needed by our program, E < 10 11 GeV. In the top part of the figure we compare the 
results of MRS R1-R2 and CTEQ 4M (with their different values of A described before) to the cross section used in 
the TIG model. In the bottom part we show the same comparison, done just with the new MRST, with its different 
A's (in all these figures cross sections increase for increasing values of A). 

All these cross sections were calculated using the MNR program, with the 'calibration' described in Sect. O, up to 
the NLO contribution. We can see in the figure that all our cross sections agree at low energies, as expected due to 
our 'calibration' at 250 GeV, and are very similar for energies up to 10 6 GeV. Beyond this energy they start showing 
their dependence on the A value and also a slight dependence on the PDF used, which was already noticed in GGV1. 
As it can be seen from both parts of the figure, the increase of the cross sections with A is evident at the highest 
energies: at the maximum energy considered the cross sections for the two extreme values of A differ by almost a 
factor of ten. 

We also notice that, for energies above 10 4 GeV, our cross sections are always considerably higher than the one 
used by TIG. As we have already explained in GGV1, TIG used an option of PYTHIA by which the gluon PDF is 
extrapolated for x < 1CP 4 with A = 0.08. In fact the TIG cross section at the highest energies shows the same slope 
of our results for A ~ 0, but is always lower than our lowest cross sections by about a factor of three. 

This can be explained only in part by the fact that the TIG cross section up to NLO is the LO result obtained 
with PYTHIA, multiplied by a constant K factor of 2, while at large energies the K factor (see GGV1 for details) 
is actually larger than 2 by about 10-15%. The bulk of the difference is however due to the different evaluations of 
the cross sections, even at LO, done by the MNR routines (our method) and directly by PYTHIA (approach used by 
TIG). 

Our results for the prompt fluxes are shown in Figs. 2-5, for MRS R1-R2, CTEQ 4M and MRST. 

In Figs. 2 and 3 we show the Ef -weighted vertical prompt fluxes Effa, calculated to NLO, for muons and muon- 
neutrinos, together with the fluxes from TIG, both from prompt and conventional sources (dotted lines). The flux 
of electron-neutrinos is practically the same as that of muon-neutrinos. Fig. 4 describes the spectral index of the 
differential fluxes, defined as ag = —d In <pg/d In Eg. 

The effects of the different extrapolations of g(x) to x < 10~ 5 (see Sect. ||) are noticeable at Eg > 10 5 GeV. In 
Figs. 2 and 3, the ^-weighted fluxes increase with A: they can differ by up to two orders of magnitude at the highest 
energy considered, 10 9 GeV, for the two extreme choices of A. This behavior is similar for all the PDF's considered. 

The A dependence of the fluxes can also affect the energy at which the prompt contribution dominates over the 
conventional sources: this is particularly true for the muon fluxes as it can be seen in Fig. 2; for the + P M fluxes this 
effect is less important (see Fig. 3) and it doesn't exist for the v e + i7 e fluxes, for which the conventional contribution 
is much lower. Apart from these differences due to the A values, charm decay dominates over conventional sources at 
Efj, ^ 10 6 GeV for muons, E v > 10 5 GeV for muon-neutrinos, and E Uc > 10 4 GeV for electron-neutrinos. 

We also see that all our fluxes for A ~ are similar to those of TIG at energies above 10 6 GeV. We have already 
mentioned that TIG used a very low value of A, A = 0.08. It is remarkable that, for these low values of A, we obtain 
similar final fluxes in spite of the differences of the two simulations and of the total cross sections already noted in 
Fig. 1. 

We can also compare our fluxes to those of the recent PRS results 0. As we have already noticed in GGV1, for 
intermediate values of A our results are very similar to the PRS ones. From Fig. 3, for example, we see that our fluxes 
for the A = 0.3 case (calculated with MRST) are close to the corresponding PRS results shown in Fig. 8 of Ref. [Q, 
calculated with CTEQ 3M and A ~ 0.3. Our results are lower than the PRS by 30 — 50% at the highest energies, 
which is probably due to the PDF's used and to the different approach of the two groups. 

Regarding the dependence of the spectral index ag on the slope A of the gluon PDF, we notice in Fig. 4 that, for 
all four PDF's, above about 10 6 GeV the differences in slope between the A = and A = 0.5 fluxes is about 0.5, 
suggesting that the spectral index is ag(Eg) = bg(Eg) — A, namely, 

MEt) ~ Ef at{El) = Ef bt(Et)+x , (14) 

where bg(Eg) is an energy dependent coefficient, that can be read off directly from the A = curve (bg(Eg) is the 
spectral index for A = 0). We will justify this result in Sect. [v| Due to this linear dependence of the spectral index 
on A, given a model which specifies the function bg(Eg), the value of A could be determined through a measurement 
of any of the <pg fluxes at two different energies. We will study in detail this possibility elsewhere jlq] . 

Here we only comment on the typical rates in a km 3 detector. It can be estimated from the curves of Fig. 2 that 
the number of prompt atmospheric muons traversing a km 3 detector from above would be over 100 per year around 
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a muon energy of 1 PeV, decreasing rapidly to less than 1 per year above 100 PeV. In this energy range there is a 
concrete possibility of detecting these prompt muons. Notice that the intensity of the prompt muon flux depends 
critically on the value of A, suggesting still another way to estimate A through the measurement of the fluxes. 

In Fig. 5 we study the dependence of the prompt fluxes on the PDF for fixed values of A. We summarize our 
previous results for A = (left) and for A = 0.5 (right), and compare them again to TIG. The figures on the top show 
the .fif -weighted fluxes, those on the bottom the spectral indices. As we already noticed in GGV1, the dependence 
on the PDF is not strong, all fluxes are very similar. This indicates that our procedure for the 'calibration' of our 
simulation with different PDF's (described in Sect. |J) is good. There are, however, some differences between the 
PDF's: in some cases (especially for A = 0) the results of MRS R2 and CTEQ 4M are very similar and higher than 
those of MRS Rl and MRST (also very close to each other). The maximum difference between all these fluxes is at 
the level of 30 to 70% at high energies. 

We want here to remark once more that our A = fluxes are very close to that of TIG at energies above I0 6 GeV 
(and also below 10 3 GeV, but the prompt fluxes are not important at these low energies). For increasing values of 
A, our results are higher than TIG, even by two orders of magnitude for A = 0.5 and at the highest energies. From 
the bottom part of the figure we notice that also the spectral indices are almost independent of the PDF used. This 
indicates that the linear dependence between ae and A of Eq. ( |l4| ) is not affected by the choice of the PDF and again 
might be used to determine the value of A. We will return on this analysis in more details in another paper fisj| . 



V. ANALYTIC INSIGHT 



In this section we first find the characteristic values of the partonic momentum fractions in the cosmic ray nucleus 
and in the nucleus in the atmosphere, and then derive the linear relation between the slope of the atmospheric muon 
(or neutrino) fluxes and the slope of the gluon parton distribution function. 

We first show that the characteristic values of the partonic momentum fractions of the incoming cosmic ray parton, 
x\, and of the target parton belonging to a nucleus in the atmosphere, x%, are respectively, 

x x ^ lO" 1 x 2 (E/10 GeV)- 1 (15) 

where E is the energy of the incoming nucleon (a proton in this paper) in the atmosphere reference frame. Precisely 
because of the small value of X2, for the relevant energies E > 10 4 GeV the gluon density g(x2) is much larger than 
the density of quarks, which we, thus, neglect in these analytic arguments. 

Let us first consider the charm flux at production d<f> c {E c , X) / 'dX , defined as the rate of c quark production^] 
per unit area, unit depth and unit charm energy (E c in the atmosphere reference frame) in the interactions of the 
attenuated nucleon flux c/>n(E,X) with the air nuclei in the atmospheric layer between X and X + dX. To obtain 
d<p c (E c , X) I dX for a layer of transverse area A and height \dh\, we simply multiply the c production rate per air 
nucleus (which equals the incoming nucleon flux at depth X times the cross section for N + A — > c + Y , where Y 
stands for "anything" and N is simply a proton p in our study) by the number of nuclei A in the layer (which is 
A\dh\riA(h)) and divide the result by the transverse area A and the layer thickness dX = p(h)\dh\. We find 

d4c{E c ,X) ^n A (h) [°° ^ da( P A^cY;E,E c ) 

~^^ = ?W/ Bc dEME > x) dE c • (16) 

We assume that the charm production cross section simply scales as A, which is expected when it is much smaller 
than the total inelastic cross section. In this case, the sum over A becomes trivial, and we have (u is the atomic mass 
unit) 

= <^*> - 1 r dE ME ^ X) MVN ^cY-^E^ 



dX u ,/ E ' dE c 



This is what we compute in our simulations (we use our 'single' mode), only the production of a c quark is calculated. Then 
the result is multiplied by two to include the contribution of the antiquark (see H for details) . 
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In these analytical considerations, we assume a simple power law for the primary flux and an energy independent 
attenuation length.^ With these approximations, the attenuated primary flux reads (see Eqs. 8-13) 

4> n {e,x) = ^{x)e-^-\ (18) 

where 4>{X) = <po exp(— X/Ajy). Substituting this approximate expression for the attenuated primary flux and chang- 
ing the integration variable from E to xe — E c /E in Eq. (p"7j), we find 

dME c ,X) _ 4>{X) x f 1 7 da(pN-*cY;x E ,E c ) 



dX u c Jo dx 

The differential cross section da(pN — > cY)/dxE is given in terms of the partonic differential cross section d&ij/dxE 
(where % and j are partons belonging to the projectile 1 and the target 2 respectively), and the PDF's fl(xi,p F ) and 
/|(ar 2 ,/i|) as 

da{pN — > cY) sr-^ f -y. 2 2 2 



dxi 



= E/ dx 1 dx 2 f!(x 1 , f ? F )ff(x 2 ,^ F )^. (20) 



Here x% and x 2 are the momentum fractions of the projectile and target partons. Mangano et al. ||] give the 
partonic cross section in terms of functions hij as 

c lFk = iP^ M^2,P,Mk,Mf), (21) 

where and -E c are the momentum and energy of the produced c quark, and, in the notation of Ref. p = 4mj?/s, 
t x = 1 - t\ — t 2 , n = (k ■ pi/pi ■ pa)) t 2 = (k ■ p 2 /p\ ■ p 2 ) and s = (pi +P2) 2 , while pi and p 2 are the projectile 
and target parton momenta respectively, p\ = x\P\^p 2 = x 2 P 2 . The hats indicate quantities in the partonic center of 
mass (those without hats are in the lab. frame at rest with the atmosphere). 

In the partonic center of mass frame, the projectile and target parton momenta are 

Pi=^,0,0,^, p 2 =^, ,0,-^, k=(E c ,0,kr,k), (22) 

and we have 

E c + k 2E C 

r 2 = 7=^, t x = 1 (23) 

V s V s 

Then, after integration over azimuthal angles, 

d 3 k d 3 k 

= — ^ = 2nd E c dk = TxsdT 2 dr x . (24) 

E c E r 



The kinematic bounds m c < E c < Vs/2 and \k\ < y E^ — m 2 c fix the integration domains of r 2 and t x . Using 
p = irric/s, we get (1 — -y/1 — p)/2 < t 2 < (1 + y/1 — p)/2 and < t x < 1 — t 2 — [p/Ar 2 ). We can use the relation 

E c k-P 2 k-p 2 

XE = -ET = -5 5" =a; l =XiT 2 , (25) 

E Pi ■ P 2 pi • p 2 
to write the differential cross section in dxE as 



2 The dependence of Ajv on E is actually very mild. In fact the whole factor e X / A «( B ) behaves like E & with (3 ~ 0.1 for 
E > 10 6 GeV and /3 even smaller for E < 10 6 GeV. Including this contribution in our analytic argument would just mean to 
replace 7 with 7 + (3 everywhere, i.e. the total spectral index would become 7 + l + /3~3.1 instead of 3.0, for energies above 
the knee at E — 5 10 6 GeV. This slight change can actually be seen in our results of Fig. 7b (see the description of that figure). 
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d&ij 
cLxe 



XlT 2 ). 



The bound X\x% — s/s > 4m 2 /2m p E — Acxe {m-p is the proton mass, m p ~ 1 GeV), where we define 



2m p E c ' 



(26) 



(27) 



implies that Xi and x 2 have a minimum lower bound larger than zero. In fact, x\ > Aex F /x 2 > 4ex^; (since x 2 < 1). 
Taking xi as the independent variable, then 4exs < xi < 1 and Acxe/xi < x% < 1< We now change the order of the 
integrations, in order to perform the integration in xe before the integrations in x\, x 2 and t 2 . 

The integration over xe in Eq. ( |l9| ) then becomes trivial, amounting to the replacement of x E by xJt^ , except for 
the necessary changes in the integration domains which become < xi,x 2 ,t 2 < 1 and < xe < (xiX 2 /e)T 2 (l — t 2 ). 
For the 8{xe — x\t 2 ) in Eq. ( p6| ) to yield a non-zero result, we need to take < x\t 2 < (xix 2 /e)T 2 (l — t 2 ), which 
means that t 2 < 1 — (e/xv) , and given that t 2 > 0, this means x 2 > e. This leads to a factorization of the xi and £2 
integrations as follows: 



cLxe x E 



da(pN -> cY) _ ^a 2 s {p R ) 



dxE 



(28) 



E 



where the functions £,j are defined as 



dr 2 t? +1 



1—V—T2 



dr x hi j(r x ,T2, 4vt 2 ,/j, R ,{i F ) 



(29) 



and the argument v is v = e/x 2 (to rewrite the integration in t 2 we noticed that p/4r 2 



The functions h 



= h^(r 2 ,p)5(r x ) 



0{a s ). We will take only gluons as partons from now on, thus 



are given by hij(T x ,T 2 , p, Rr, Vf] 
fl(x,v F ) = ff(x,(i F )=g(x,(i F ). 

The function using h gg at the Born level, is shown in Fig. 6a for 7 = 1.7 and 2 (corresponding to the spectral 
indices 7 + 1 of the primary flux above and below the knee) . In the same figure we see that the maximum of £ ffg (v) is 
at v ~ 0.1, namely x 2 ~ 10 e. However, given that g(x 2 ,p^,) is a sharply increasing function with decreasing x 2 (i.e. 
for increasing v at fixed E c ), the maximum of the product g(x 2 , p 2 F )t,gg{ v ) is always to the right of the maximum of 
C, gg {v), at v > 0.1. Therefore, the integral in x 2 in Eq. (28) is dominated by the values of x 2 of order e, namely 



x 2 



GeV 



Returning to Eq. (28), the integral in x\ shows that large values of x\ will be dominant since xjg(xi) 



(30) 

7-A-1 

for small x, where the exponent is positive, since 7 = 1.7 or 2, while < A < 0.5 (thus 7 — A — 1 > 0). To see 
more precisely what range of x\ dominates the integral, it is necessary to prove two statements. The first is that 
t 2 = xe/xi < 1, due to kinematical constrains, therefore x\ > x F - The second is that the characteristic value of Xe 
is 0.1, namely that the c-quark is mainly produced with 0.1 of the proton energy 



E c = O(0.1 E), 



(31) 



With respect to the kinematical limit on t 2 , as we already mentioned, t 2 = xe/x\ < 1 — v, and we obtained as a 
kinematical constraint that e < v = e/ x 2 < 1 (since x 2 goes from e to 1). Thus, t 2 < 1 — e < 1, since e is always larger 
than zero. Another way of obtaining this bound is the following. Since the partonic processes involved are gg — * cc 
or gg — > ccg, then y/I > 2(E c ) max and due to m c ^ 0, (/c) ma x < (E c )m«, therefore r 2 < 2(£ c ) max /v / I < 1. 

That in fact E c = O(0.1 E) is clearly demonstrated in Fig. 6b, which shows the function x E [da j dxE) normalized by 
the total c-production cross section. Thus we have proven that the dominant range of x x in Eq. (28) is X\ > O(0. IE) 



and also, combining together Eq. (|3C|) and Eq. (|3l|), our statement in Eq. ( |l5[) about x 2 . 

Even if we have not yet included gluon shadowing in our calculations, we want to point out that this effect might 
only be important for the target gluon (given that x 2 is very small) but it is not important for the gluons in the 
projectile (given that x\ ^ 0.1). This means that the uncertainties on the composition of cosmic rays will not affect 
the results through shadowing effects. 
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As a summary of our arguments we can say that, due to the incoming flux being rapidly falling with increasing 
energy of the primary, only the charm quarks produced with a large fraction of the incoming energy, E c ~ 0.1-B, 
count in the charm flux at production, and those highly energetic c quarks come from projectile partons carrying a 
large fraction of the incoming momentum x\ > xe — 0.1. On the other hand, because typical partonic center of mass 
energies V§ are close to the cc threshold, 2m c — 2 GeV (since the cross section decreases steeply with increasing y/l), 
while the total center of mass energy squared is s = 2m p E (with m p the proton mass, m p ~ 1 GeV), the product 
Xl x 2 = s/s = 4m 2 J(2m p E) ~ GeV/E. This shows that x 2 ^ (GeV/Exx) ~ GeV/0.1B. 

We now derive the dependence on A of the muon and neutrino fluxes for a simple power law primary flux. 

We can explain first the dependence on A of the spectral index of d<p c /dX at large energies, and then, using this 
result, the dependence on A of the spectral indices of atmospheric muons and neutrinos. To start with, we notice 
that the integral in Eq. (28) depends on the charm energy E c only through the presence of the parameter e in the 
integration on X2- To approximately perform this integration at large energies, let us replace g{ X 2) — %2 ~~ m 
Eq. (28) and take £(e/xz) ~ £max (namely develop £ in powers of v = e/x2 and keep only the constant term) then 

dx 2 g(x 2 ) C( 3- ] — Cmax / dx 2 Xz^ 1 ■ (32) 




Since t«l, this integral is well approximated by £ max e~ /A, for all A ^ 0. Better approximations to the function 
£ give similar results. For example, approximating the function £ by two power laws, one above and another below 
the maximum, which is at about X2 = 5e (£ = Cmax(^2/5e) 21 for X2 between e and 5e and £ = Cmax(5e/a;2) for X2 
between 5e and 1), the integral in Eq.(p2]) becomes Cmax(5e) _A /(0.9 + 1.7A — A 2 ). Thus the essential dependence of 
e~ A is maintained. Recalling that e = m~/(2 m p E c ), Eq. (|l9| ) is proportional to E x , and the same is true for Eq. 
(|32"|), therefore 

^E e ,X)~Eri-™ . (33) 

The charm production function d<j) c (E c , X)/dX, calculated numerically, is shown in Fig. 7a for a typical X — 
57.12 g/cm 2 (h = 20 km). We are using here the PDF MRS Rl with the three related values of A = 0, A(R1), 0.5. 
We clearly see here that the slope at E c > 10 5 GeV depends on the extrapolation of the gluon PDF at x < 10~ 5 . This 
is one order of magnitude lower in energy than in Fig. 1 for the total cross section. This reflects the fact mentioned 
above that the characteristic charm energy is E c = 0(0. IS). Fig. 7b shows that, as predicted analytically, the slopes 
(the negative of the spectral index in our notation) of the charm fluxes at production depend almost linearly on A. In 
fact, in Fig. 7b, we can see that the logarithmic slopes of the A = and A = 0.5 fluxes differ precisely by 0.5, above 5 
10 6 GeV (namely, above the knee) to about 10 9 GeV (the maximum energy at which our fluxes are reliable, given that 
we take 10 11 GeV as the maximum incoming proton energy E). In fact, the slope of the A = flux in that interval is 
about -3.1 to -3.2, while that of the A = 0.5 flux is about -2.6 to -2.7. Above the knee, the primary spectrum goes as 
E s with d ~ (-7 - 1 - 0.1) = -3.1, where we have also included the 0.1 contribution coming from the independence 
of An (see footnote in previous discussion), thus the charm spectrum, (in the energy range 10 7 GeV < E c < 10 9 GeV) 
goes approximately as E^ +x as expected from Eq. ((33|). 

Using the definition of the leptonic fluxes in terms of the charm spectrum at production d(f) c /dX, we can now find 
the dependence of the spectral index of muon and neutrino fluxes with A. For example, the differential flux M of 
muons with energy E p (/i stands here for /i + or /i~) is 



poo poo 

«/v(£ M ) = 2 dX dE, 

J Xr\ J En 



d(j> c (E c ,X) 



dE^ 



(34) 



(ipfj. has, thus, units of [1/ cm s sr GeV]). Here the factor of 2 accounts for the muons produced by c and the last 
square bracket is the number of muons of energy E p produced at sea level by the cascades, each cascade initiated by 
a c quark of energy E c at a depth X . 

Our results above indicate that we can write the atmospheric charm spectrum at production as (see Eq.(|33|)) 
d(f>c(E c ,X) jdX ~ F{X)E^ 1 ^ 1JrX with F(X) a function independent of energy. Replacing this form for d(f> c (E c , X) jdX 
in Eq. ( p4| ) and multiplying and dividing by -E^ 7 ~ 1+A w e can write </> M as 



MEJ = 2E^- 1+X / dXF(X) / dE c (-f- 

JXn JE U 



1 1+A r dN^c-twEcE^X) 
dE„ 



(35) 
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We can argue that in so far as the values of the parent charm quark energy E c and the daughter lepton energy 
are not very different, the dependence of the integral on A (and on 7) should be mild. In this case, from Eq. (|35|), we 
find that the spectral index of the muon (and similarly of the neutrino) flux contains A as a term, i.e. 

<M25„) * f{E^,\)E^- 1+x = i5-MBM.7,A)+A , (36) 

where the dependence of the functions f{E u , 7, A) and 6^(^,7, A) on A and 7 should be mild. This justifies the 
results shown in Figs. 4 and 5, presented in Sect. IV, showing all the spectral indices obtained using all our PDF's. 

Finally we examine the deviations from linearity of the relation between the spectral index a. a and the gluon PDF 
slope A. In Fig. 8a we show directly the relation between A and a^, using the values coming from our simulation 
for the MRST case already presented in Fig. 4, but now plotting them for fixed energy E^. We show two examples, 
for E u = 1 PeV 1 10 PeV , where our points indicate a good agreement with the linear relation between at and A of 
Eq. Q). 

The mild dependence on A of the functions bi(\) = ae + A can be seen in Fig. 8b, where we show the percentage 
difference [b e (X) - b e (0)]/be(0) for the different values of A = 0-0.5 with the MRST PDF. It is evident that, in the 
range where our theoretical arguments are applicable (for £„ > 10 6 GeV) the bi{X) functions differ only by 2 — 3% 
for different A values, namely they are almost independent of A, given one particular PDF. This analysis confirms 
the validity of Eq. ( |14| ) , which leads to the possibility of obtaining information on A at small parton fractions x not 
reachable in experiments, through the measurement of the fluxes. We will study this possibility in more detail in a 
future paper UH]. 



VI. CONCLUSIONS 



The actual next-to-leading order perturbative QCD calculations of charm production cross sections, together with 
a full simulation of the atmospheric cascades, were used to obtain the vertical prompt fluxes of neutrinos and muons. 

We have analyzed the dependence of the atmospheric fluxes on the extrapolation of the gluon PDF at very low x, 
which is related to the value of the parameter A. This was done using four different sets of PDF's: MRS Rl, MRS 
R2, CTEQ 4M and MRST, with variable A in the range 0-0.5. 

The charm production cross sections and the final lepton fluxes depend critically on A for leptonic energies Ei > 
10 5 GeV, which correspond to x < 10 -5 GeV. We found that the fluxes vary up to almost two orders of magnitude 
at the highest energy considered, 10 9 GeV, for the different A's in the allowed interval; on the contrary, for fixed A, 
the results don't depend much on the choice of the PDF. 

For the lowest values of A (A ~ — 0.1) our fluxes are very close to those of TIG ||, confirming that the very low 
flux prediction is mostly due to a low value of A (A TIG ~ 0.08). For higher values of A (A ~ 0.2 — 0.5) our results 
are in the bulk of previous predictions and, in particular, for A ~ 0.3 they are very close to a recent semi- analytical 
calculation j?J done with a similar value of A. 

We have also considered the dependence of the spectral index of the final fluxes on the parameters of the model. 
From both, computer simulations and analytical considerations, we find that the spectral index at of atmospheric 
leptonic fluxes depends linearly on A as in Eq. (|l4|). 

This suggests the possibility of obtaining bounds on A in "neutrino telescopes" for small values of x not reachable 
in colliders, if the spectral index of leptonic atmospheric fluxes could be determined by these telescopes. We will 
investigate this possibility in detail in the future Jl8[] . 
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FIGURE CAPTIONS 

Fig. 1 Total cross section for charm production a C c, up to NLO, for our different PDF's and A values, compared to 
that used by TIG §]. Top panel: MRS R1-R2 and CTEQ 4M; bottom panel: MRST (cross sections increase 
with A). 

Fig. 2 Prompt muons: i? 3 -weighted vertical fluxes at NLO, compared to the TIG || conventional and prompt fluxes 
(dotted lines). We show results using the four PDF's MRS Rl, MRS R2, CTEQ 4M and MRST. 

Fig. 3 Prompt muon-neutrinos: i? 3 -weighted vertical fluxes at NLO, compared to the TIG || conventional and 
prompt fluxes (dotted lines). We show results using the four PDF's MRS Rl, MRS R2, CTEQ 4M and MRST. 

Fig. 4 Prompt muons: spectral index of the NLO vertical fluxes for the four PDF's MRS Rl, MRS R2, CTEQ 4M 
and MRST. 

Fig. 5 Dependence of prompt fluxes and their spectral index on the PDF at fixed A: left side A = 0, right side A = 0.5. 

Fig. 6 (a) The function Q g (v) at the Born level for 7 = 0, 1.7 (below the knee) and 7 = 2 (above the knee), (b) 
Flux- weighted charm production spectra Xg-j^p a * severa l beam energies (using MRS Rl, A(R1)). 

Fig. 7 (a) NLO charm production function E^d(j) c (E cl X)/dX (PDF MRS Rl); (b) its spectral index 
—d\n[d(f) c (E c ,X)/dX}/d\nE c . These results are for a height h — 20 km, corresponding to a vertical depth 
X = 57.12 g/cm 2 (similar results are obtained for other heights). 

Fig. 8 (a) Relation between the slope A of the gluon PDF and the muon spectral index a M at fixed muon energy, (b) 
Non-linearities in this relation. Here bi(X) = ai(\) + A and we use the MRST PDF. 
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